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ABSTRACT 

Econometricians  must  choose  between  many  methods  for  estimating  p,  the 
autocorrelation  coefficient,  in  a  first  order  autoregressive  process.  This  thesis  examines 
the  performance  of  four  estimators  in  a  Monte  Carlo  simulation.  The  methods 
examined  are  Durbin- Watson,  Beach-MacKinnon,  Theil-Nagar  and  Prais-Winsten. 
The  autocorrelation  coeficient,  p,  was  varied  from  .2  to  .9  and  each  method  provided 
estimates  of  p  and  P,  the  regression  coefficient,  for  1000  replications.  The  results 
presented  here  are  similar  to  those  found  in  previous  comparisons.  Specifically, 
Ordinary  Least  Squares  was  found  to  be  an  efficient  estimator  of  P  when 
autocorrelation  is  present  only  to  a  slight  degree.  Of  the  four  estimators  examined,  the 
performance  of  Theil-Nagar  proved  superior  in  estimating  both  p  and  p  for  small 
values  of  the  correlation  coeficient.  Beach-MacKinnon,  on  the  other  hand,  while 
containing  a  large  bias  in  the  estimation  of  p,  is  the  more  efficient  estimator  of  P  for 
large  values  of  p. 


TABLE  OF  CONTENTS 


I.  INTRODUCTION 8 

A.  BACKGROUND    8 

B.  PROBLEM  STATEMENT    8 

C       ESTIMATORS 9 

D.  SIMULATION 9 

E.  MEASURE  OF  EFFECTIVENESS   9 

II.  ESTIMATION .  10 

A.  GENERAL 10 

B.  PROCESS  10 

C.  METHODS  OF  ESTIMATION   12 

1.  Generalized  Least  Squares  Estimation 12 

2.  Estimated  Generalized  Least  Squares 12 

3.  Maximum  Likelihood  Estimation 14 

III.  COMPARISON 16 

A.  GENERAL 16 

B.  PREVIOUS  COMPARISONS 16 

C.  MODEL  AND  DATA  GENERATION 17 

D.  VALIDATION 17 

E.  SIMULATION 17 

F.  MEASURES  OF  EFFECTIVENESS  18 

IV.  RESULTS  AND  CONCLUSIONS 20 

A.  GENER.\L 20 

B.  SAMPLE  SIZE  20   20 

C.  SAMPLE  SIZE  50   22 

D.  SUMMARY   22 

APPENDIX  :    PROGRAM  LISTINGS  27 

LIST  OF  REFERENCES 47 


INITIAL  DISTRIBUTION  LIST , 48 


LIST  OF  TABLES 

L  ESTIMATES  OF  AUTOCORRELATION  COEFFICIENT     ; 18 

II.  ESTIMATES  OF  AUTOCORRELATION  COEFFICIENT     21 

III.  DATA  PRESENTED  IN  FIGURE  4.1    23 

IV.  EFFICIENCY  OF  REGRESSION  COEFFICIENT 

ESTIMATES    24 

V.  DATA  PRESENTED  IN  FIGURE  4.2   25 

VI.  EFFICIENCY  OF  REGRESSION  COEFFICIENT 

ESTIMATES    26 


LIST  OF  FIGURES 

4.1  Estimated  mean  square  error  of  p  vs.  p  (sample  size  =  20) 23 

4.2  Estimated  mean  square  error  ofp  vs.  p  (sample  size  =  50) 25 


I.  INTRODUCTION 

A.  BACKGROUND 

Autocorrelation  exists  in  a  regression  model  when  the  error  terms  are  no  longer 
independent  but  are  correlated.  In  the  examination  of  time  scries  data  autocorrelation 
is  a  common  phenomenon  and  can  lead  to  problems  if  Ordinary  Least  Squares  (OLS) 
estimation  procedures  are  used.  The  purpose  of  this  thesis  is  to  examine  and  compare 
four  different  estimates  of  the  autocorrelation  coeflTicicnt,  p,  the  estimation  of  which  is 
essential  to  the  resolution  of  OLS  deficiencies.  The  four  estimators  to  be  examined  are 
the  Durbin- Watson,  Theil-Nagar,  Beach-MacKinnon,  and  Prais-Winsten. 

B.  PROBLEM  STATEMENT 

In  the  standard  regression  model  y=XP  +e,  y  is  a  Txl  vector  of  observations  of 
a  dependent  variable,  X  is  a  TxK  design  matrix  and  P  is  a  Kxl  vector.  The  variable  e 
is  a  Txl  vector  of  unobservable  random  errors  with  E(o)=0  and  covariancc  matrix, 
E{ee')=cr  Ij.  Thus,  in  the  standard  model,  the  random  vector  e  contains  elements 
which  are  pairwise  uncorrclated  with  identical  means  and  variances.  In  the  presence  of 
autocorrelation  this  strong  assumption  is  violated.  That  is,  the  error  terms  arc  no 
longer  independent  but  arc  correlated.   The  regression  model  becomes, 

y^  =  X^P  +  ej  t=l,2,....,T  (eqn  I.l) 

where      e^.  =  p  e^-.j  +  v^, 
E(v^)  =  0,  and 
E(vv')  =  (T^I  . 

This  is  known  as  a  first  order  autoregressive  or  AR(1)  process.    As  illustrated  by 

equation  1.1,  e^.  is  expressed  linearly  in  terms  of  the  e^j    and  another  random  error 

term  v..    The  assumption  of  zero  mean  and  constant  variance  provides  v.  with  all  the 

nice  properties  of  e^.  in  the  standard  model.    This  process  may  occur  for  a  variety  of 

reasons,  some  of  which  are: 

/  Omitted  expUmatorv  variables.  If  a  correlated  cxplanatorv  variable  has  been 
excluded  from  the  clesign  matrix  its  exclusion  will  be  rcllect'ed  in  the  correlation 
of  the  random  variable  e. 

2  Mispecification  of  the  mathematical  Jonn  of  the  model.  If  the  wrong  mathematical 
relationship  is  chosen  the  values  ol  e  may  be  dependent. 


3  Interpolations  in  the  statistical  observations.    If  the  observational  data  is  smoothed 
autocorrelation  may  result. 

4  Mispecification  of  the  true  random  error.    Dependence  among  the  error  terms  may 
occur  naturally.   [Ref.  l:p.   204j 

Utilizing  OLS  to  estimate  the  regression  coelTicient,  P,  in  the  presence  of  an  AR(1) 
process  can  lead  to  problems.  Generally,  there  are  two  consequences  to  consider.  The 
first  is  that  the  OLS  estimator  of  the  coefficients  will  be  unbiased  but  will  not  be  very 
efficient.  The  second  consequence  is  that  the  OLS  variance  estimator  is  biased.  For 
these  reasons  it  is  useful  to  investigate  other  methods  to  estimate  P  [Ref  2:p.   439]. 

C.  ESTIMATORS 

When  p  is  known,  the  process  is  easily  accounted  for  using  Generalized  Least 
Squares  or  Weighted  Least  Squares  methods  [Ref  3].  However,  the  usual  situation  is 
that  p  is  unknown  and  must  be  estimated.  A  number  of  methods  have  been  proposed 
to  estimate  p  and  properly  account  for  OLS  deficiencies  in  estimating  p.  Chapter  2  will 
develop  the  four  estimators  mentioned  above  and  examine  the  autocorrelation  process. 

D.  SIMULATION 

Each  of  the  estimators  considered  here  have  the  same  asymptotic  properties 
therefore  any  decision  on  which  one  to  use  must  be  based  on  small  sample  analysis  and 
Monte-Carlo  evidence.  Therefore,  a  simulation  will  be  created  in  which  the  data  is 
generated  according  to  guidelines  presented  in  previous  studies  with  equation  1.1  as  the 
model.  The  actual  values  of  p  will  be  varied  from  .2  to  .9.  The  four  estimation 
techniques  will  then  provide  estimates  of  p  and  p  for  1000  replications. 

E.  MEASURE  OF  EFFECTIVENESS 

To  provide  an  indication  of  which  estimator  performs  best  the  mean  square  error 
of  both  p  and  P  will  be  estimated  for  each  estimator.  Prior  results  for  different  sets  of 
estimators  indicate  that  no  one  estimator  will  prove  superior  over  the  entire  range  of  p 
but  that  one  or  two  may  out  perform  the  others  over  specific  intervals. 


II.  ESTIMATION 

A.  GENERAL 

This  chapter  attempts  to  explore  the  theory  behind  both  the  first  order  process 
and  four  estimators  developed  to  properly  account  for  it.  Three  of  these 
(Durbin-Watson,  Thiel-Nagar,  and  Prais-Winsten)  arc  categorized  as  estimated 
generalized  least  squares  estimators.  The  fourth  (Beach-MacKinnon)  is  a  maximum 
likelihood  estimator. 

B.  PROCESS 

The  first  order  process  can  be  written  as 

y^  =  X^P  +  e^  t=   1,2,...,T  (eqn2.1) 

where  e^.  =  p  e^._j  +  v^, 

E(v^)  =  0, 

E(vv')  =  G^l, 

E(v^2)=<yy^,  and 

E(VjV^)  =  0  for  S3:t  . 

The  parameter  p  is  generally  unknown  and  along  with  P  must  be  estimated.  The 
Statistical  properties  of  the  random  error,  v,  listed  in  equation  2.1  are  identical  to  those 
listed  for  e  in  the  general  linear  model.  The  statistical  properties  of  e  under  these  new 
assumptions  are  quite  different.   Judge  [Ref  4:p.   438]  shows  that 

E(ej)  =  2  p'E{v^.i)  =  0  (eqn  2.2) 

and 

E(e^^)=(j^^=<y^^/(\-p)^.  (eqn  2.3) 

The  covariance  between  errors  s  periods  apart  is  no  longer  zero  and  is  given  by 

E(eje^.3)=E(ej^^e^)  =  (p^(T\)/(I-p2)  .  (eqn  2.4) 
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The  covariance  matrix  for  e  is  now  easily  written  as 


a>  =  E(ee')  = 


(eqn  2.5) 


<tV(i-p^) 


1         P 
P         1 


T-1     T-2 


or  utilizing  the  following  convention, 


,T-1 
.T-21 


0)  =  g/'V 


(eqn  2.6) 


where  T  = 


1/(1-P^) 


2 
P        P 


,T-1 
,T-2 

»T-3 


T-1     T-2     T-3 
pi   1  pi  z,  pi  J 


Thus,  the  assumptions  made  about  the  error  term,  e,  in  the  standard  linear  model 
no  longer  hold  for  the  autoregressive  case.  Specifically,  due  to  autocorrelation  the  error 
covariance  matrix  is  no  longer  written  as  a  I  but  is  now  (T     T  . 


When  an  attempt  is  made  to  perform  a  least  squares  fit  to  the  data  in  the  presence 

of  an  AR(1)  process  there  are  two  problems  to  consider. 

1    The  least  squares  estimator  P  =  (X'X)    X'y  will  be  unbiased  but  will  not  be  veiy 
efilcicnt. 


2    The  least  squares  covariance  matrix  d  (X'Xj      with  a 
will  be  a  biased  estimator  of  the  variance  of ' 


(y-Xb)'(y-Xb)/(T-K) 


In  the  presence  of  positive  autocorrelation  Judge  [Ref  4:p.    439]  notes  that  with 


>N 


OLS  estimation  the  bias  of  the  standard  error  of  P  will  very  likely  appear  as  an 
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underestimate.  Park  and  Mitchell  [Ref.  5:p.  16]  warn  that  OLS  seriously 
underestimates  the  variance  of  P  for  p  >  0.4.  This  understatement  makes  the 
estimates  themselves  appear  much  more  significant  than  they  actually  are  and  makes 
hypothesis  testing  of  the  slope  coefficients  unreliable. 

C.       METHODS  OF  ESTIMATION 

I.    Generalized  Least  Squares  Estimation 

When  apriori  information  is  available  about  T,  the  most  convenient  estimate 
.         A    . 
for  the  regression  coefficient,   P,  is   obtained  by  applying  least   squares  estimation 

techniques  to  the  transformed  model, 


Y*  =  X-p  +  e* 
where  Y*  =  PY 
X-=  =  PX 
e*  =  Pe. 


(eqn  2.7) 


The  transformation  matrix  P  is  the  TxT  matrix 


P  = 


-p 

0 

0 


0  0 

1  0 

-p  1 

0  -p 


0      0 


where  PT  =  T^^ 


-P 


This  method  is  known  as  the  Generalized  Least  Squares  (GLS)  estimation. 
2.    Estimated  Generalized  Least  Squares 

The  usual  case  is  that  p  is  unknown  and  must  be  estimated.  Once  an  estimate 
for  p  (p)  is  computed  one  can  substitute  p  into  the  P  matrix  and  proceed  with  the  GLS 
method  outlined  above.  This  is  known  as  Estimated  Generalized  Least  Squares 
(EGLS)  estimation.  The  computational  form  of  the  alternative  estimators  for  p 
discussed  are  as  follows: 
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a.   Durbin-  Watson 
The  statistic 

d  =  j;  (e^-e^.i)2/i;^^2       ^=  ^^^^  (,qn  2.8) 

where  e^  =  y^-Xj.)) 


is  often  used  to  test  for  first  order  autoregressive  errors.  As  the  number  of 
observations  (T)  increases  it  can  be  demonstrated  that  d  approaches  the  least  squares 
estimator  of  p  or 

^=   1-  (d/2)  .  (eqn  2.9) 

The  Durbin- Watson  statistic  is  provided  by  most  least  squares  computer  packages  and 
is  very  easy  to  use.  It  also  is  an  example  of  a  two-stage  estimator.  That  is,  it  first 
estimates  the  correlation  parameter  and  then  uses  this  estimate  to  compute  the 
generalized  least  squares  estimates  for  p. 

b.  Theil-Nagar 

A  modification  of  the   Durbin-Watson  estimator  suggested  by  Henri 
Theil  and  A.  L.  Nagar  is 

p-  =  (T^(I-(d/2))  +  K^)  /(T^  -  K^).  (eqn  2.10) 

Theil  and  Nagar  claim  that  this  estimator  is  an  improvement  over  Durbin-Watson  if 
the  first  and  second  differences  of  the  explanatory  variables  are  small  when  compared 
to  their  corresponding  ranges  [Ref  6].  Like  Durbin-Watson,  it  also  is  a  two-stage 
estimator. 

c.  Prais-  IV  ins  ten 

A  minimum  sum  of  squares  approach  to  estimating  p  yields, 

.  T  T-, 

P  =  I   Vt-l/SV       t=l,...,T  (eqn  2.11) 


where  e^  =  yt'^t^ 
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This  estimator  can  be  employed  in  both  a  two  step  and  an  iterative  procedure.    This 
paper,  however,  considers  only  the  following  iterative  form: 

1.  Set'p  =  0. 

2.  Transform  the   variables   in   accordance   with   the   transformation  matrix   and 

equation  2.7. 

„  A 

3.  Calculate  the  least  squares  estimate  of  p  conditional  on  p. 

4.  Calculate  the  estimate  of  p  conditional  on  P  by  using  equation  2.11. 

5.  If  the  absolute  difference  in  p  from  the  previous  iteration  is  suificiently  small  (less 
than  0.00001)  stop.    If  not  go  to  step  2.   [Ref  7:p.   2J 


3.   Maximum  Likelihood  Estimation 

A  maximum  likeUhood  (ML)  estimator  is  the  value  of  6  which  maximizes  the 

value  of  the  likelihood  function  L(0).    Under  the  assumption  that  Y  has  a  multivariate 

■J 
normal  distribution  with  mean  XP  and  covariance  matrix  <7  T,  the  likelihood  function 

is 

L(P,p,(t2)=  C-  (l/2(T\)(y-XP)'T-^(y-XP)  (cqn  2.12) 

where  C  =  -{li2)\n(p-^  +  (l/2)ln  (l-p^)  . 


The  ML  estimators  for  P,  p,  and  (Ty    are  those  values  for  which, 

aL/5P  =  0,    aL/ap  =  0,    dLlda^^  =  Q   .  (eqn2.13) 


Solutions  to  equations  2.13  are  very  difficult  to  derive.  Beach  and  MacKinnon 
[Ref  8:p.  54]  use  an  ML  estimator  for  <3^  and  substitute  into  equation  2.12.  The 
result  is  the  concentrated  likelihood  function, 

L(P.p)=  K-(T/2)ln((y-Xp)'T-ky-XP)(l-p2)l/'r)  (eqn  2.14) 

where  K  =  (T/2)ln(T)-(T/2)  . 

They  suggest  maximizing  L(p,p)  with  respect  to  P  with  p  held  constant  and  then  to 
maximize  with  respect  to  p  with  P  held  constant.  An  algorithm  to  derive  this  ML 
estimate  is 
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1.  Set  p  =  0. 

2.  Transform  the  variables  in  accordance  with  equation  2.7. 

3.  Calculate  the  least  squares  estimate  of  P  conditional  on  p. 

4.  Calculate  the  ML  estimate  of  p  conditional  on  P  by  solving  a  cubic  equation  of 
the  untransformed  residuals,  (see  [Ref  8]  for  details) 

5.  If  the  absolute  difference  in'p  from  the  previous  iteration  is  sufficiently  small  (less 
than  O.OOOOI)  stop.  If  not,  go  to  step  2  [Ref  7].  (Note:  The  same  procedure  was 
employed  for  iterative  Prais-Winsten  method  except  that  equation  2. 1 1  was  used  to 
estimate  p.) 

This  is  not  a  comprehensive  listing  of  all  available  estimators  for  a  first  order 
process.   Other  estimators  are  listed  in  Judge  [Ref  4]. 
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III.  COMPARISON 

A.  GENERAL 

The  finite  sampling  properties  of  the  estimators  Hsted  here  have  not  been  derived. 
Choice  of  which  estimator  to  use  might  be  based  on  evidence  obtained  from  Monte 
Carlo  simulations.  This  chapter  explains  a  simulation  used  and  provides  a  synopsis  of 
comparisons  reported  in  the  literature. 

B.  PREVIOUS  COMPARISONS 

There  have  been  a  number  of  studies  of  estimators  for  p  .  Each  has  concluded 
that  OLS  has  serious  deficiencies  in  the  presence  of  autocorrelation.  The  majority  of 
these  papers  have  settled  on  two  points.  First,  particularly  in  small  sample  sizes 
(T<50)  it  is  best  to  use  estimators  that  consider  all  T  observations.  Rao  and 
Grilitches  concluded  that  using  estimators  such  as  Cochrane-Orcutt  that  ignore  the 
first  observation  can  lead  to  a  substantial  loss  of  efficiency  [Ref  9:p.  269J.  These 
results  were  further  substantiated  by  Beach  and  MacKinnon.  In  an  attempt  to  develop 
a  computationally  efficient  algorithm  to  maximize  the  likelihood  function  they 
discovered  (for  p=  0.6,  0.8,  0.99)  significant  gains  in  efficiency  to  be  made  by 
employing  the  first  observation.  Some  of  these  gains  are  in  the  neighborhood  of  700 
percent  [Ref  8;p.  55].  Park  and  Mitchell  concluded  that  retention  of  this  first 
observation  substantially  reduces  the  risk  of  collincarity  as  p  approaches  0.9  [Ref  5:p. 
10].  Kobayashi  verified  theoretically  the  experimental  results  of  Park  and  Mitchell.  By 
computing  the  asymptotic  variances  of  several  estimators  he  demonstrated  that  the  loss 
of  efficiency  of  the  Cochrane-Orcutt  method  was  due  primarily  to  ignoring  the  first 
observation.   [Ref  10:p.  951]. 

The  second  point  is  that  the  Prais-Winsten  solution  techniques  outperform  many 
comparable  estimators  of  the  correlation  parameter.  Spitzer  concluded  that 
Prais-Winsten  "appeared  to  be  the  best  of  all  the  two  stage  estimators."  [Ref  11  :p. 
44].  Park  and  Mitchell  in  a  later  study  comparing  Beach-MacKinnon  with  the  iterative 
Prais-Winsten  estimator  concluded  that  the  iterative  Prais-Winsten  performs 
"appreciably  better  in  estimating  the  autocorrelation  coefficient  p"  [Ref  7:p.   5]. 
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Although  there  were  no  studies  found  specifically  comparing  the  four  estimators 
presented  here,  each  has  demonstrated  a  superiority  to  OLS  in  the  presence  of  a  first 
order  process. 

C.  MODEL  AND  DATA  GENERATION 

Equation  2.1  was  utilized  as  the  model  with  the  first  term  in  the  vector  e 
generated  in  the  following  fashion, 

ei=  Vi/(l-p)l/2  (eqn3.1) 

In  order  to  conform  with  previous  comparisons,  the  data  utilized  in  this 
experiment  is  identical  to  that  used  in  Beach  and  MacKinnon  [Rcf  8].  Two  sample 
sizes  of  20  and  50  observations  were  used.  The  untrended  explanatory  variable,  X,  was 
drawn  from  N(0,  0.0625)  and  the  random  error,  v^,  was  drawn  from  N(0,  0.0036). 
Although  autocorrelation  in  theory  may  be  positive  or  negative,  in  econometric  data  it 
is  almost  always  positive  [Rcf  l:p.  201].    For  this  reason  p  was  varied  from  0.2  to  0.9. 

D.  VALIDATION 

The  data  generation  program  was  checked  to  ensure  the  normality  of  c  using  the 
Chi  Square  Goodness  of  Fit  test.  The  normality  assumption  was  accepted  at  a  0.3684 
level.  Finally,  in  order  to  ensure  each  estimator  performed  properly  the  random 
portion  of  the  model,  specifically  the  random  variable  V,  was  removed.  This  allowed 
the  estimators  to  function  in  a  deterministic  fashion.  Data  were  then  generated  and 
submitted  to  each  estimator  for  values  of  p  equal  0.2,  0.6,  0.8.  The  results  are 
presented  in  Table  I  ,  illustrate  that  the  estimators  are  functioning  properly. 

E.  SIMULATION 

For  each  run  the  values  of  the  regression  coefficients,  Pq  and  Pi,  were  set  to  1 
and  1.  The  variables  X^.  and  v^  were  drawn  from  the  normal  distributions  discussed 
earlier.  The  dependent  variable  y^  was  calculated  using  equation  2.1.  Since  the 
ultimate  objective  was  to  generate  residuals  to  send  to  the  four  estimation  routines,  a 
regression  was  then  performed  of  y  on  X  and  residuals  calculated  using, 

"^^  =  yfX^p^  t  =  1,2,...,T.  (eqn3.2) 
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TABLE  I 

ESTIMATHS  OF  AUTOCORRELATION  COEFFICIENT 

P 

DW 

TN 

PW 

BM 

.2 

.19 

.19 

.19 

.19 

.6 

.59 

.60 

.58 

.60 

.8 

.78 

.80 

.77 

.80 

The  values  of  the  residuals  were  then  sent  to  each  estimation  routine.  Estimates  of 
P  and  p  were  determined  for  values  of  p  equal  to  .2,  .3,  .4,  .5,  .6,  .7,  .8,  and  .9.  Each 
estimate  was  replicated  1000  times  for  the  sample  sizes  of  20  and  50. 

F.       MEASURES  OF  EFFECTIVENESS 

In  order  to  compare  the  performances  of  the  estimators,  two  MOE's  were  used. 
The  mean  square  error  (MSE)  of  p  was  estimated  for  each  estimator.  '1  his  represents 
the  expected  squared  error  made  in  estimating  p.  The  following  computational  form  of 
MSE  was  used, 

V(P-Pi)/1000       i=  1,...,1000.  (eqn  3.3) 

/■-/ 


A 

The  successive  values  of  MSI:  of  p  were  then  plotted  against  the  actual  p  to 
examine  performance  over  the  range  of  p. 

The  second  MOE  examined  the  relative  efficiencies  of  the  regression  coefficient  as 
defined  in  [Ref  7:p.    7].   A  ratio  of  MSE  of  p  for  a  particular  estimate  to  the  MSE  of  P 
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for  the  OLS  estimate  allows  the  examination  of  the  relative  gains  in  using  particular 
techniques  over  OLS.  Since  the  proper  estimation  of  {5  is  paramount  the  efficiency  of  p 
is  predetermined  to  be  the  most  important  MOE. 
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IV.  RESULTS  AND  CONCLUSIONS 

A.  GENERAL 

The  major  emphasis  of  this  thesis  was  to  examine  the  performance  of  four 
estimators  of  the  autocorrelation  coefficient,  p,  for  a  first  order  autoregressive  process. 
The  estimators  examined  were  Durbin-Watson,  Theil-Nagar,  Prais-Winsten,  and 
Beach-MacKinnon. 

A  Monte-Carlo  simulation  was  performed  for  the  following  values  of  p:  .2,  .3,  .4, 
.5,  .6,  .7,  .8,  and  .9.  Each  run  was  replicated  1000  times  for  sample  sizes  of  20  and  50. 
The  results  are  recorded  in  Table  II  .  Irrespective  of  sample  size,  each  of  the  methods 
underestimate  the  true  value  of  p  but  as  the  number  of  observations  is  increased  from 
20  to  50  the  bias  reduces.  As  was  expected,  no  one  estimator  uniformly  outperforms 
the  others.  In  both  sample  sizes,  the  two  stage  estimators  (Durbin-Watson  and 
Theil-Nagar)  achieve  better  results  for  small  p.  As  the  value  of  p  increases,  the 
iterative  methods  (Prais-Winsten  and  Beach-MacKinnon)  perform  best.  With  T=20 
this  transition  occurs  at  p=  .6  while  at  50  observations  it  occurs  earlier  at  p=  .A. 

The  discussion  of  the  results  will  be  divided  into  two  sections.  The  measures  of 
effectiveness,  as  defined  in  Chapter  3  will  first  be  applied  to  the  simulation  results  for 
T  =  20.  This  will  be  followed  by  an  identical  approach  when  the  sample  size  is 
increased  to  50. 

B.  SAMPLE  SIZE  20 

Since  performance  of  an  estimator  is  roughly  indicated  by  its  mean  and  variance, 

A 

mean  square  error  (MSC)  of  each  p  over  the  entire  sample  size  was  estimated.  'I he 
results  of  these  calculations  arc  presented  in  Table  III  along  with  a  plot  of  MSC  of  p 
versus  actual  values  of  p  in  Figure  4.1  .  They  again  indicate  that  the  Theil-Nagar  and 
Durbin-Watson  estimators  arc  better  for  smaller  values  of  p  (p<.6)  and  as  p  increases 
the  Prais-Winsten  p  emerges  as  the  best.  On  the  basis  of  Figure  4.1  alone, 
Beach-MacKinnon's  performance  is  clearly  inferior.  However,  in  examining  the 
eniciency  of  each  estimator  in  Table  IV,  Beach-MacKinnon  proves  to  be  the  most 
efficient  in  estimating  P  over  the  widest  range  of  p.  The  tie  in  Figure  4.1  between 
Theil-Nagar  and  Durbin-Watson  is  resolved  in  Table  IV  with  Theil-Nagar  proving  to 
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TABLE  II 

ESTIMATES  OF  AUTOCORRELATION  COEFITCIENT 

Sample  Size 

20 

P 

DW 

TN 

PW 

BM 

.2 

.158 

.162 

.113 

.107 

.3 

.234 

.239 

.205 

.193 

.4 

.310 

.316 

.296 

.279 

.5 

.385 

.392 

.387 

.365 

.6 

.460 

.467 

.478 

.450 

.7 

.533 

.542 

.567 

.535 

.8 

.603 

.613 

.655 

.617 

.9 

.667 

.681 

.741 

.697 

Sample  Size  50 

P 

DW 

FN 

PW 

BM 

.2 

.173 

.161 

.170 

.178 

.3 

.279 

.253 

.270 

.270 

.4 

.360 

.340 

.360 

.359 

.5 

.451 

.432 

.463 

.453 

.6 

.544 

.523 

.559 

.547 

.7 

.630 

.610 

.651 

.640 

.8 

.726 

.700 

.747 

.731 

.9 

.812 

.796 

.839 

.820 
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be  uniformly  more  efficient  than  Durbin-Watson.    Table  IV  also  demonstrates  that  for 
p  =  .2  OLS  is  at  least  as  efTicient  as  three  of  the  four  estimators. 

C.  SAMPLE  SIZE  50 

The  results  of  the  MSE  calculations  for  T=50  are  recorded  in  Table  V  along 
with  a  plot  of  MSE  of  p  versus  the  actual  values  of  p  in  Figure  4.2.  The 
Durbin-Watson  and  Theil-Nagar  estimators  again  perform  the  best  for  smaller  values 
of  p  (p  <  .4)  and  as  p  increases  the  Beach-McKinnon  and  Prais-Winsten  estimators  of 
p  contain  the  smallest  MSE. 

Once  again  even  though  the  Prais-Winsten  p  has  a  smaller  MSE  than 
Beach-MacKinnon,  Table  VI  illustrates  that  Bcach-McKinnon  is  a  uniformly  more 
efficient  estimator  of  the  slope  coefficient.  For  the  smaller  values  of  p  (p  <  .4) 
Theil-Nagar  is  more  efficient  than  Durbin-Watson.  Table  VI  also  illustrates  that  OLS 
is  at  least  as  efficient  as  any  of  the  other  estimators  when  p  is  small. 

D.  SUMMARY 

As  was  found  in  previous  studies  when  autocorrelation  is  present  only  to  a  slight 
degree  (p<.2)  the  OLS  estimator  provides  an  efficient  estimate  for  the  regression 
coefficient,  p.  As  the  process  becomes  more  significant  however,  all  the  estimators 
outperform  the  OLS  solution.  In  both  sample  sizes  the  performance  of  Theil-Nagar 
and  Durbin-Watson  are  nearly  identical  with  respect  to  the  MSE  of  p.  However,  when 
efficiency  of  the  slope  coefficient  estimate  is  examined,  Theil-Nagar  proves  to  be  the 
better  2  stage  estimator.  Park  and  Mitchell  [Ref  7:p.  4J  found  that  Prais-Winsten 
performs  better  in  estimating  p.  The  results  presented  here  tend  to  dispute  that 
finding.  For  while  Prais-Winsten  has  a  uniformly  smaller  MSE  of  p, 
Beach-MacKinnon  provides  the  most  efficient  estimator  of  P  .  Spitzer,  on  the  other 
hand  [Ref  ll:p.  44],  which  ranked  two  stage  estimators  as  being  the  best  for  values  of 
p  between  .2  and  .5,  mirrors  the  results  produced  here.  Apriori  knowledge  of  the 
neighborhood  of  p  will  be  helpful  in  selecting  the  appropriate  estimation  method.  For 
both  sample  sizes  Theil-Nagar  appears  to  be  the  best  for  small  values  of  p. 
Beach-MacKinnon,  while  containing  a  larger  bias  for  p  than  does  Prais-Winsten,  is  a 
much  more  efficient  estimator  of  the  slope  coefficient  for  larger  values  of  p. 
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Figure  4.1     Estimated  mean  square  error  ofp  vs.  p  (sample  size  =  20). 


TABLE  III 

DATA  PRnSlINTCD  IN  FIGURl:  4.1 

Sample  Si 

/e  20 

P 

MSEDW 

MSETN 

MSEPW 

MSEBM 

.2 

.7494 

.7485 

.8557 

.8594 

.3 

.6268 

.6244 

.7004 

.7115 

.4 

.5156 

.5124 

.5621 

.5787 

.5 

.4159 

.4126 

.4400 

.4604 

.6 

.3278 

.3250 

.3342 

.3566 

.7 

.2519 

.2495 

.2433 

.2662 

.8 

.1891 

.1860 

.1698 

.1916 

.9 

.1407 

.1357 

.1141 

.1331 
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TABLE  IV 

rrriciiiNCY  or  regrcssiox  (:()i:rFicinNT  estimates 

SaiTi 

pie  Si/e  20 

p 

MSEP  (DW) 

MSEP  (TN) 

MSEP  (PW) 

MSEP(»M) 

MSEP(()I,S) 

MSEp  (OLS) 

MSEP  (OLS) 

MSEP  (OLS) 

.2 

1.004 

.9794 

1.035 

1.041 

.3 

.9228 

.8967 

.9442 

.9515 

.4 

.8218 

.7929 

.8325 

.8342 

.5 

.7082 

.6751 

.7024 

.6959 

.6 

.5864 

.5484 

.5652 

.5515 

.7 

.4610 

.4207 

.4329 

.4135 

.8 

.3359 

.3020 

.3093 

.2870 

.9 

.2253 

.2087 

.2077 

.1892 
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Figure  4.2     Estimated  mean  square  error  ofp  vs.  p  (sample  size  =  50). 


TABLE  V 

DATA  PRESHN  ri:D  IN  FIGURE  4.2 

Sample  Si 

ze50 

P 

MSEDW 

MSETN 

MSEPW 

MSEBM 

.2 

.7010 

.6766 

.7055 

.7065 

.3 

.5500 

.5399 

.5653 

.5578 

.4 

.4383 

.4196 

.4396 

.5578 

.5 

.3298 

.3151 

.3056 

.3156 

.6 

.2358 

.2262 

.2116 

.2215 

.7 

.1578 

.1526 

.1357 

.1449 

.8 

.0906 

.0942 

.0773 

.0851 

.9 

.0500 

.0509 

.0360 

.0417 
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TARI.F  VI 

HFFIClIiNCY  OF  RFGRFSSION  COFFI  ICIIINT  FSTIMATFS 

Sampl 

e  Size  50 

p 

MSFP{DW) 

MSFP(TN) 

MSFP(PW) 

MSFP  (DM) 

MSFP  (OLS) 

MSFp(OFS) 

MSFP  (OLS) 

MSFP  (OLS) 

.2 

1.073 

1.041 

1.046 

1.058 

.3 

.9985 

.9482 

.9714 

.9562 

A 

.8850 

.8255 

.8635 

.8255 

.5 

.7452 

.6859 

.7282 

.6825 

.6 

.5920 

.5420 

.5870 

.5406 

.7 

.4366 

.4020 

.4453 

.4020 

.8 

.2889 

.2690 

.3067 

.2700 

.9 

.1589 

.1505 

.1738 

.1505 
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APPENDIX 
PROGRAM  LISTINGS 

This  appendix  contains  listings  of  the  programs  utilized  in  the  analysis  performed 
herein.  All  of  the  functions  are  written  in  FORTR.'\N  and  contain  the  necessary 
documentation.  The  Monte  Carlo  simulation  was  performed  using  the  Advanced 
Simulation  and  Statistics  Package  [Ref  12]  developed  by  P.  A.  Lewis.  Since  the 
package  only  allows  for  the  simultaneous  comparision  of  3  estimators,  2  functions  were 
developed  for  each  sample  size.  The  fust,  SIMS  generates  estimates  for 
Durbin-Watson,  Theil-Nagar,  and  Prais-Winstcn  for  a  sample  size  of  20.  SI  MSA 
meanwhile,  generates  estimates  for  Beach-MacKinnon  for  the  identical  sample  size. 
Routines  for  Durbin-Watson  and  Theil-Nagar  were  included  in  SI  MSA  to  ensure  the 
results  were  comparable  to  SIMS.  SI  MSB  and  SIMSC  perform  in  a  similar  fashion  for 
sample  size  of  50  and  therefore  were  not  included.  The  Advanced  Simulation  and 
Statistics  Package  computes  the  mean  square  error  of  'p'  for  each  estimator 
automatically.  The  mean  square  error  for  the  p  estimates  was  accomplished  by  the 
MSEB  function. 


SIMS 
DIMENSION  EHAT(20) 
COMMON  /MYDATA/  K,T,ANS,Y1,X 
COMMON  /DATAl/  IX1A,RH0 
REAL*4  Y(5000),YMIN,YMAX,PMEAN(3) 
CHARACTER*80  T1,T2,T3 

INTEGER  N, M, NE(8),L,D,RG,SEI,SVS, NEST, NSR, 1X1,1X2, 1X3 
EXTERNAL  DATGEN,  DURWAT,  BEAMAC,  PRAWIN,  LSEB,  DCALC,  TRANSF 
EXTERNAL  LNORM,SIMTBD,GMPRD 
NR=20 
T=20 
K=2 
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OPEN  (UNIT=19,FILE='M0NICA') 
C    OPEN  (UNIT=21,FILE='MARGE') 
C    OPEN  (UNIT=51,FILE='AMBR0SE') 
C    OPEN  (UNIT=41,FILE='DAT2') 

OPEN  (UNIT=51,FILE='DAT3') 

READ  (19,*)  ANS 
10    READ(19,*,END=999)  N,M,L,D,RG,SEI ,SVS,NEST,NSR 

READ(19/)YMIN,YMAX 

READ(19,*)  (NE(I),I=1,L) 

READ(19,120)  1X1,1X2,1X3 
120   F0RMAT(I5,1X,I5,1X,I5) 

READ  (19,115)  Tl 
115   FORMAT(A80) 

READ(19,115)  T2 
■   READ  (19,115)T3 

READ(19,*)  (PMEAN(I),I=1,3) 

READ(19,*)  RHO 

READ(19,61)IX1A 
61    F0RMAT(I5) 
C 

C     CALL  FOR  SIMTBD 
C 

CALL  SIMTBD  (IXl , 1X2, IX3,Y,N,M,NE, L,D,NSR,RG,SEI ,SVS, 
*YMIN,YMAX,NEST,DATGEN,DURWAT,T1,DATGEN,BEAMAC,T2,DATGEN,PRAWIN,T3, 
*PMEAN) 

GO  TO  10 
999   WRITE(6,*)'END  OF  DATA  INPUT' 

STOP 

END 

r     ********)!<****  DATA  GENERATION  SUBROUTINE  *********y<)i(***yr*5ir****)it 

p      ***************************************  y(*)'c)«e)ic*yc*ycyc***)'c)<(yc*!»()'t*yc)'c**)>t* 

SUBROUTINE  DATGEN  (1X1 ,EHAT,NR) 
DIMENSION  BHAT(2) ,YSTAR(20) ,R2(20) ,U(20) , 
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*E(20),YHAT(20),EHAT(20),XSTAR(20,2) 
*,Y1(20),X(20,2),V(20) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

COMMON  /DATAl/  IXIA.RHO 

INTEGER  IX1,IX1A,NR 
C 
C 

C    GENERATE  THE  RANDOM  ERROR 
C 

CALL  SNOR  (IX1,U,NR,1,0) 
C 
C    ADJUST  THE  VARIANCE  OF  R.  E.  lAW  BEACH  AND  MACKINN0N(1978) 

DO  38  1=1, T 

V(I)=U(I)*.06 
38    CONTINUE 
C 

C    GENERATE  THE  ERROR  FOR  THE  STAND  LINEAR  MODEL 
C 

E(1)=V(1)/(1-(RHO**2))**0.  5 

DO  31  J=2,T 

E(J)=RH0*E(J-1)+V(J) 

31  CONTINUE 
C 

c 

C     GENERATE  THE  EXPLANATORY  VARIABLES  lAW  RAO  AND  GRILITCHES  (1969) 
C 

DO  32  1=1,20 
X(I,1)=1 

32  CONTINUE 

C     CHANGE  IXl  IN  ORDER  TO  AVIOD  COLLINEARITY 
C     IX1A=IX1+19 

CALL  SNOR(IX1A,R2,NR,1,0) 

DO  33  J=l,20 

X(J,2)=R2(J)*.25 

33  CONTINUE 

29 


c 
c 

C    THE  TRUE  BETA  EQUALS  1,1 

C 

C 

C     GENERATE  THE  INDEPENDENT  VARIABLE 

C 

DO  35  1=1,20 

Y1(I)=(X(I,1)+X(I,2))+E(I) 
35    CONTINUE 
C 

C    GENERATE  THE  LEAST  SQUARES  ESTIMATOR  FOR  BETA 
C 

CALL  LSEB(X,Y1,BHAT) 
C     PRINT  LSEB  TO  A  FILE 

IF(ANS  .EQ.  2)  WRITE(61,201)BHAT 
201   F0RMAT(F11.8,2X,F11.8) 
C 

C    GENERATE  YHAT 
C 

DO  100  1=1,20 

YHAT(I)=X(I,1)*BHAT(1)+X(I,2)*BHAT(2) 
100   CONTINUE 
C 
C 

C    GENERATE  EHAT 
C 

DO  50  1=1,20 

EHAT(I)=Y1(I)-YHAT(I) 
50    CONTINUE 
C 
C 

RETURN 

END 
r  ******************yf**  DURBIN  WATSON  ************************ 
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C    THIS  FUNCTION  COMPUTES  THE  DURBIN-WATSON  ESTIMATE  OF  RHO 

REAL  FUNCTION  DURWAT  (EHAT,NR,WI) 

DIMENSION  EHAT(20),X(20,2),Y1(20),XSTAR1(20,2),YSTAR1(20),BHAT1(2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

CALL  DCALC  (EHATJ.D) 

DURWAT=l-D/2 
C 

CALL  TRANSF(X,Y1, DURWAT, XSTAR1,YSTAR1) 

CALL  LSEB  (XSTAR1,YSTAR1,BHAT1) 

IF  (ANS  .EQ.  1  )  WRITE(21,701)  BHATl 
701   F0RMAT(F11.8,2X,F11.8) 
C 
C 

c 

RETURN 

END 
C 
C 

p      ************>!(*******   TucTi  MAPAR   ****************************  A* 

C    THIS  FUNCTION  COMPUTES  THE  THEIL-NAGAR  ESTIMATE  OF  RHO 
C     REAL  FUNCTION  THENAG  (EHAT,NR,WI) 

DIMENSION  EHAT(20),YSTAR2(20),XSTAR2(20,2),BHAT2(2) 
*,Y1(20),X(20,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

CALL  DCALC  (EHAT,T,D) 

THENAG=((T**2)*(l-D/2)+K**2)/(T**2-K**2) 

CALL  TRANSF(X,Y1, THENAG, XSTAR2,YSTAR2) 

CALL  LSEB  (XSTAR2,YSTAR2,BHAT2) 

IF  (ANS  .EQ.  1  )  WRITE(31,801)  BHAT2 
801   F0RMAT(F11.8,2X,F11.8) 

RETURN 

END 

p     ***********************   PRATS  WTNSTFN   ************************* 

C    THIS  FUNCTION  COMPUTES  THE  PRAIS-WINSTEN  ESTIMATE  OF  RHO 
REAL  FUNCTION  PRAWIN(EHAT,NR,WI) 
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DIMENSION  EHAT3(20) ,YHAT3(20) ,YSTAR3(20) ,BHAT3(2) , 
*EHAT(20),XSTAR3(20,2) 
*,Y1(20),X(20,2) 
COMMON  /MYDATA/  K,T,ANS,Y1,X 
N=0 

RH03=0 
98    N=N+1 

CALL  TRANSF  (X,Y1 ,RH03,XSTAR3,YSTAR3) 
CALL  LSEB  (XSTAR3,YSTAR3,BHAT3) 
C     GENERATE  YHAT3 
DO  83  1=1,20 

YHAT3(I)=X(I,1)*BHAT3(1)+X(I,2)*BHAT3(2) 
83    CONTINUE 

DO  4  1=1, T 

EHAT3(I)=Y1(I)-YHAT3(I) 

4  CONTINUE 
C 

RHONUM=0 
RHODEN=0 
DO  5  1=2, T 

RHONUM=RHONUM+( EHAT3( I )*EHAT3( I-l) ) 

5  CONTINUE 
C 

DO  6  1=2, T-1 

RH0DEN=RH0DEN+(EHAT3(I)**2) 

6  CONTINUE 

PRAWIN=RHONUM/RHODEN 
C    CHECK  FOR  PRAWIN  WHICH  ARE  OUT  OF  BOUNDS 
IF(PRAWIN.GE.  1)THEN 

PRAWIN=0. 99999 
ELSE  IF  (PRAWIN.  LE.-1)THEN 

PRAWIN=-0.  99999 
END  IF 
C     COMPARISION  OF  RH03  AND  PRAWIN  IF  DIFF  .  LT.  0.0001  THEN  END 
IF(ABS(RH03-PRAWIN).  GT. .  0001)THEN 
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RH03=PRAWIN 
GO  TO  98 
ELSE 

PRAWIN=PRAWIN 
END  IF 
C     IF  (ANS  .EQ.  1  )  WRITE(41,901)  BHAT3 
COl   F0RMAT(F11.8,2X,F11.8) 
RETURN 
END 
C 

C    THE  FOLLOWING  SUBROUTINES  AID  IN  THE  COMPUTATION  OF  THE  FOUR 
C     ESTIMATORS  OF  RHO. 

p      ytAyc**********   SUBROUTINE  LSEB   *****************yt**)k*ycyf******* 

C     SUBROUTINE  LSEB  WILL  COMPUTE  THE  LSE  OF  B 
C 

SUBROUTINE  LSEB(X,Y1 ,BHAT) 

DIMENSION  BHAT(2),Y1(20),X(20,2),XTRNSP(2,20),XI(2,2),H(2,20), 
*XPRIX(2,2) 
C     X  TRANSPOSE 

DO  40  1=1,20 

DO  41  J=l,2 

XTRNSP(J,I)=X(I,J) 
41      CONTINUE 
40    CONTINUE 
C    MULTIPLY  X  TRANSPOSE  AND  X 

CALL  GMPRD(XTRNSP,X,XPRIX,2,20,2) 
C     CALCULATE  INVERSE  OF  X  PRIME  X 

DETR=1/(XPRIX(1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 

XI(1,1)=DETR*XPRIX(2,2) 

XI(1,2)=DETR*(-XPRIX(1,2)) 

XI(2,1)=DETR*(-XPRIX(2,1)) 

XI(2,2)=DETR*XPRIX(1,1) 
C    MULTIPLY  INVERSE  AND  TRANSPOSE 

CALL  GMPRD(XI,XTRNSP,H,2,2,20) 

DO  99  1=1,2 
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BHAT(I)=H(I,1)*Y1(1)+H(I,2)*Y1(2)+H(I,3)*Y1(3) 
*+H(I,4)*Yl(4)+H(I,5)*Yl(5) 

*+H(I,6)*Yl(6)+H(I,7)*Yl(7)+H(I,8)*Yl(8)+H(I,9)*Yl(9) 
*+H(I,10)*Yl(10)+ 

*H(I,11)*Y1(11)+H(I,12)*Y1(12)+H(I,13)*Y1(13)+H(I,14)*Y1(14) 
*+H(I,15)*Yl(15)+ 

*H(I,15)*Y1(16)+H(I,17)*Y1(17)+H(I,18)*Y1(18)+H(I,19)*Y1(19)+ 
*H(I,20)*Y1(20) 
99    CONTINUE 

RETURN 

END 

C     SUBROUTINE  DCALC  WILL  COMPUTE  THE  DURBIN  STATISTIC  D 
C 

SUBROUTINE  DCALC(EHATJ,D) 

DIMENSION  D1(20),D2(20),EHAT(20) 

DNUM=0 

DDEN=0 

DO  1  1=2, T 

D1(I-1)=(EHAT(I)-EHAT(I-1))**2 
DNUM=DNUM+D1(I-1) 

1  CONTINUE 

DO  2  J=1,T 

D2(J)=EHAT(J)**2 
DDEN=DDEN+D2(J) 

2  CONTINUE 

D=DNUM/DDEN 

RETURN 

END 
C 
p     *****************  SUBROUTINE  TRANSF  ********************** 

C 

C     SUBROUTINE  TRANSF  IS  DESIGNED  TO  TRANSFORM  THE  X'S  AND  Y'S 
C    ACCORDING  TO  THE  LEAST  SQUARES  RULE. 

SUBROUTINE  TRANSF(X,Y1,RH0HAT,XSTAR,YSTAR) 
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DIMENSION  Y1(20),YSTAR(20),X(20,2),XSTAR(20,2) 

K=2 

T=20 
C    Y  TRANSFORM 

YSTAR(1)=((1-(RH0HAT**2))**0. 5)*Y1(1) 

DO  7  1=2,20 

YSTAR(I)=Y1(I)-(RH0HAT*Y1(I-1)) 
7     CONTINUE 
C     X  TRANSFORM 

XSTAR(1,1)=(1-(RH0HAT**2))**0.  5 

DO  9  J=2,K 

XSTAR(1,J)=((1-(RH0HAT**2))**0.5)*X(1,J) 
9     CONTINUE 

DO  11  L=2,T 

XSTAR(L,1)=1-RH0HAT 

11  CONTINUE 

DO  12  1=2, T 
"DO  13  J=2,K 

XSTAR(I,J)=X(I,J)-RH0HAT*X(I-1,J) 
13      CONTINUE 

12  CONTINUE 

RETURN 
END 
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SIMSA 
C    THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  RUN  COMPUTE  THE  FOLLOWING 
C     ESTIMATORS  (DW  TN  BM)  FOR  A  SAMPLE  SIZE  OF  20 

DIMENSION  EHAT(20) 

COMMON  /MYDATA/  KJ,ANS,Y1,X 

COMMON  /DATAl/  IX1A,RH0 

REAL*4  Y(5000),YMIN,YMAX,PMEAN(3) 

CHARACTER*80  T1,T2,T3 

INTEGER  N, M, NE(8),L,D,RG,SEI,SVS, NEST, NSR, 1X1,1X2, 1X3 

EXTERNAL  DATGEN,  DURWAT,  THENAG,  BEAMAC,  LSEB,  DCALC,  TRANSF 

EXTERNAL  LNORM,SIMTBD,GMPRD 

NR=20 

T=20 

K=2 
C 
C 

c 

OPEN  (UNIT=19,FILE='M0NICA') 

OPEN  (UNIT=51,FILE=' AMBROSE') 

READ  (19,*)  ANS 
10    READ(19,*,END=999)  N,M,L,D,RG,SEI ,SVS, NEST, NSR 

READ(19,*)YMIN,YMAX 

READ(19,*)  (NE(I),I=1,L) 

WRITE  (22,105)  (NE(I) , 1=1 , L) 
105   F0RMAT(8I4) 

READ(19,120)  1X1,1X2,1X3 
120   F0RMAT(I5,1X,I5,1X,I5) 

READ  (19,115)  Tl 
115   F0RMAT(A8a) 

READ(19,115)  T2 

READ  (19,115)T3 

READ(19,*)  (PMEAN(I),I=1,3) 

READ(19,*)  RHO 

READ(19,61)IX1A 
61    F0RMAT(I5) 
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c 

C     CALL  FOR  SIMTBD 
C 

CALL  SIMTBD  (IXl ,IX2,IX3,Y,N,M,NE,L,D,NSR,RG,SEI ,SVS, 
*YMIN,YMAX,NEST,DATGEN,DURWAT,T1,DATGEN,THENAG,T2,DATGEN,BEAMAC,T3, 
*PMEAN) 
GO  TO  10 
999   WRITE(6,*)'END  OF  DATA  INPUT' 
STOP 
END 

p     *************  DATA  GENERATION  SUBROUTINE  ********************** 
p     ***************************************************************** 

SUBROUTINE  DATGEN  (IX1,EHAT,NR) 

DIMENSION  BHAT(2),YSTAR(20),R2(20),U(20), 
*E(20),YHAT(20),EHAT(20),XSTAR(20,2) 
*,Y1(20),X(20,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

COMMON  /DATAl/  IXIA.RHO 

INTEGER  IX1,IX1A,NR 
C  _        _ 

c 

C    GENERATE  THE  RANDOM  ERROR 
C 

CALL  SNOR  (IX1,U,NR,1,0) 
C 

c 

C    GENERATE  THE  ERROR  FOR  THE  STAND  LINEAR  MODEL 
C 

E(1)=U(1)/(1-(RHO**2))**0. 5 

DO  31  J=2,20 

E(J)=RH0*E(J-1)+U(J) 
31    CONTINUE 
C 

c 
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C    GENERATE  THE  EXPLANATORY  VARIABLES  lAW  RAO  AND  GRILITCHES  (1969) 
C 

DO  32  1=1,20 
X(I,1)=1 

32  CONTINUE 

C    CHANGE  IXl  IN  ORDER  TO  AVIOD  COLLINEARITY 
C     IX1A=IX1+19 

CALL  SN0R(IX1A,R2,NR,1,0) 

DO  33  J=l,20 

X(J,2)=R2(J)*.  25 

33  CONTINUE 
C 

c 

C    THE  TRUE  BETA  EQUALS  1,1 

C 

C 

C     GENERATE  THE  INDEPENDENT  VARIABLE 

C 

DO  35  1=1,20 

Y1(I)=(X(I,1)+X(I,2))+E(I) 
35  -  CONTINUE  - 
C 

c 

C     GENERATE  YHAT 
C 

CALL  LSEB(X,Y1,BHAT) 
C    BHAT(1)=1.3 
C    BHAT(2)=1.  1 

DO  100  1=1,20 

YHAT(I)=X(I,1)*BHAT(1)+X(I,2)*BHAT(2) 
100   CONTINUE 
C 

c 

C    GENERATE  EHAT 
C 
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DO  50  1=1,20 

EHAT(I)=Y1(I)-YHAT(I) 

50    CONTINUE 

C 

C 

RETURN 

END 
r  *********************  DURBIN  WATSON  *************************** 

C 

REAL  FUNCTION  DURWAT  (EHAT,NR,WI) 

DIMENSION  EHAT(20),X(20,2),Y1(20),XSTAR1(20,2),YSTAR1(20),BHAT1(2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

CALL  DCALC  (EHAT,T,D) 

DURWAT=l-D/2 

CALL  TRANSF(X,Y1, DURWAT, XSTAR1,YSTAR1) 

CALL  LSEB  (XSTARl ,YSTAR1 ,BHAT1) 
C 
C 
C 

RETURN 

END 
C 
C 
r  ********************  THEIL  NAGAR  ****************************** 

C 

REAL  FUNCTION  THENAG  (EHAT,NR,WI) 

DIMENSION  EHAT(20),YSTAR2(20),XSTAR2(20,2),BHAT2(2) 
*,Y1(20),X(20,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

CALL  DCALC  (EHAT,T,D) 

THENAG=((T**2)*(l-D/2)+K**2)/(T**2-K**2) 

RETURN 

END 
C 
Q     ***********************  RFACH  MACKINNON  ************************* 
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REAL  FUNCTION  BEAMAC(EHAT,NR,WI) 

DIMENSION  EHAT4(20),YHAT4(20),YSTAR4(20),BHAT4(2), 
*Y1(20),EHAT(20),X(20,2),XSTAR4(20,2) 

COMMON  /MYDATA/  KJ,ANS,Y1,X 

N=0 

RH04=0 
98    N=N+1 

CALL  TRANSF  (X,Y1 ,RH04,XSTAR4,YSTAR4) 

CALL  LSEB  (XSTAR4,YSTAR4,BHAT4)  . 
C    BHAT4(1)=1.0 
C    BHAT4(2)=1.0 
C    GENERATE  YHAT4 

DO  83  1=1,20 

YHAT4(I)=X(I,1)*BHAT4(1)+X(I,2)*BHAT4(2) 
83    CONTINUE 

DO  4  1=1, T 

EHAT4(I)=Y1(I)-YHAT4(I) 
4    CONTINUE 

SUM3=0 

SUM2=0 

SUM1=0 

DO  71  1=2, T 

SUM1=SUM1+(EHAT4(I)*EHAT4(I-1)) 

71  CONTINUE 
C 

DO  72  1=2, T 

SUM2=SUM2+(EHAT4(I-1)**2) 

72  CONTINUE 
C 

DO  73  1=2, T 

SUM3=SUM3+(EHAT4(I)**2) 

73  CONTINUE 
C 

DEN0M=(T-1)*(SUM2-(EHAT4(1)**2)) 
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A=(-(T-2)*SUM1)/DEN0M 

B=(((T-1)*(EHAT4(1)**2))-(T*SUM2)-SUM3)/DEN0M 

C=(T*SUM1)/DEN0M 

SMALP=B-((A**2)/3) 

SMALQ=C-((A*B)/3)+((2*(A**3))/27) 

THETA=AC0S((SMALQ*(27**.  5))/(2*SMALP*((-SMALP)**0. 5))) 


C 

C     BEAiMAC  IS  THE  ITERATIVE  RHO  FOR  THIS  PROCEEDURE 

BEAMAC=(-2*((-SMALP/3)**0.  5))*C0S((THETA/3)+(3. 1412/3))-(A/3) 
C     CHECK  FOR  BEAMAC  WHICH  ARE  OUT  OF  BOUNDS 
IF(BEAMAC.GE.  1)THEN 

BEAMAC=0.  99999 
ELSE  IF  (BEAMAC.  LE.-1)THEN 

BEAMAC=-0.  99999 
END  IF  _        _ 

C     COMPARISION  OF  RH04  AND  BEAMAC  IF  DIFF  .  LT.  0.0001  THEN  END 
IF(ABS(RH04-BEAMAC). GT. . 0001)THEN 
RH04=BEAMAC 
GO  TO  98 
ELSE 

BEAMAC=BEAMAC 
.  END  IF 

IF  (ANS  .EQ.  2)  WRITE  (51,901)  BEAMAC 
901   FORMAT( F15.il) 
RETURN 
END 
C 

C    THE  FOLLOWING  SUBROUTINES  AID  IN  THE  COMPUTATION  OF  THE  FOUR 
C     ESTIMATORS  OF  RHO. 
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r  *************  SUBROUTINE  LSEB  ******************************* 

C     SUBROUTINE  LSEB  WILL  COMPUTE  THE  LSE  OF  B 

C 

SUBROUTINE  LSEB(X,Y1,BHAT) 

DIMENSION  BHAT(2),Y1(20),X(20,2),XTRNSP(2,20),XI(2,2),H(2,20), 
*XPRIX(2,2) 
C     X  TRANSPOSE 

DO  40  1=1,20 
DO  41  J=l,2 

XTRNSP(J,I)=X(I,J) 
41      CONTINUE 
40    CONTINUE 
C     MULTIPLY  X  TRANSPOSE  AND  X 

CALL  GMPRD(XTRNSP,X,XPRIX,2,20,2) 
C     CALCULATE  INVERSE  OF  X  PRIME  X 

DETR=1/(XPRIX(1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 
XI(1,1)=DETR*XPRIX(2,2) 
XI(1,2)=DETR*(-XPRIX(1,2)) 
XI(2,1)=DETR'^(-XPRIX(2,1)) 
XI(2,2)=DETR*XPRIX(1,1) 
C     MULTIPLY  INVERSE  AND  TRANSPOSE 
CALL  GMPRD(XI,XTRNSP,H,2,2,20) 
DO  99  1=1,2 

BHAT(I)=H(I,1)*Y1(1)+H(I,2)*Y1(2)+H(I,3)*Y1(3) 
*+H(I,4)*Yl(4)+H(I,5)*Yl(5) 

*+H(I,6)*Yl(6)+H(I,7)'^Yl(7)+H(I,8)*Yl(8)+H(I,9)*Yl(9) 
*+H(I,10)*Yl(10)+ 

*H(I,11)*Y1(11)+H(I,12)*Y1(12)+H(I,13)*Y1(13)+H(I,14)*Y1(14) 
*+H(I,15)*Yl(15)+ 

*H(I,16)*Y1(16)+H(I,17)*Y1(17)+H(I,18)*Y1(18)+H(I,19)*Y1(19)+ 
*H(I,20)*Y1(20) 
99    CONTINUE 
RETURN 
END 
p     *************  SUBROUTINE  DCALC  ***************************** 
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C     SUBROUTINE  DCALC  WILL  COMPUTE  THE  DURBIN  STATISTIC  D 
C 

SUBROUTINE  DCALC(EHAT,T,D) 

DIMENSION  D1(20),D2(20),EHAT(20) 

DNUM=0 

DDEN=0 

DO  1  1=2, T 

D1(I-1)=(EHAT(I)-EHAT(I-1))**2 
DNUM=DNUM+D1(I-1) 

1  CONTINUE 

DO  2  J=1,T 

D2(J)=EHAT(J)**2 
DDEN=DDEN+D2(J) 

2  CONTINUE 

D=DNUM/DDEN 

RETURN 

END 
C 
r  *****************  SUBROUTINE  TRANSF  ********************** 

C 

C     SUBROUTINE  TRANSF  IS  DESIGNED  TO  TRANSFORM  THE  X'S  AND  Y'S 

C    ACCORDING  TO  THE  LEAST  SQUARES  RULE. 

SUBROUTINE  TRANSF(X,Y1 ,RHOHAT,XSTAR,YSTAR) 

DIMENSION  Y1(20),YSTAR(20),X(20,2),XSTAR(20,2) 

K=2 

T=20 
C    Y  TRANSFORM 

YSTAR(1)=((1-(RHOHAT**2))**0.5)*Y1(1) 

DO  7  1=2,20 

YSTAR(I)=Y1(I)-(RH0HAT*Y1(I-1)) 
7    CONTINUE 
C     X  TRANSFORM 

XSTAR(1,1)=(1-(RH0HAT**2))**0.  5 

DO  9  J=2,K 

XSTAR(1,J)=((1-(RH0HAT**2))**0.  5)*X(1,J) 
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9     CONTINUE 

DO  11  L=2,T 

XSTAR(L,1)=1-RH0HAT 

11  CONTINUE 

DO  12  1=2, T 
DO  13  J=2,K 

XSTAR(I,J)=X(I,J)-RH0HAT*X(I-1,J) 
13      CONTINUE 

12  CONTINUE 

RETURN 
END 
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MSEB 
C    THIS  PROGRAM  IS  DESIGNED  TO  CALCULATE  THE  MEAN  SQUARE  ERROR  OF 
C    THE  BETA  VECTOR 

DIMENSION  Bl(5000), 82(5000), 83(5000), B4( 5000), B5( 5000), 86(5000), 
*B7( 5000), B8( 5000), B9( 5000), B10( 5000), BX( 5000), BY(  5000) 

OPEN  (UNIT=21,FILE='DAT1') 

OPEN  (UNIT=31,FILE='DAT2') 

OPEN  (UNIT=41,FILE='DAT3') 

OPEN  (UNIT=51,FILE='DAT4') 

OPEN  (UNIT=61,FILE='DAT5') 
C 

COUNT=1000 

READ(21,900)(B1(I),B2(I),  1=1,1000) 

CALL  MSEBET  (Bl ,B2, COUNT, XMSEDW) 

READ(31,900)(B3(I),B4(I),  1=1,1000) 

CALL  MSEBET  (B3, 84, COUNT, XMSETN) 

READ(41,900)(B5(I),B6(I),  1=1,1000) 

CALL  MSEBET  (B5,B6,C0UNT,XMSEPW) 

READ(51,900)(B7(I),B8(I),  1=1,1000) 

CALL  MSEBET  (B7, 88, COUNT, XMSEBM) 

READ(61,900)(B9(I),B10(I),  1=1,1000) 

CALL  MSEBET  (89, 810, COUNT, XMSEOLS) 
900   FORMAT  (Fll.  8,2X,F11.  8) 

WRITE(6,*)'MSEDW' 

WRITE(6,*)XMSEDW 
C 

WRITE(6,*)'MSETN' 

WRITE(6,*)XMSETN 
C 

WRITE(6,*)'MSEPW' 

WRITE(6,*)XMSEPW 
C 

WRITE(6,*)'MSE8M' 

WRITE(6,*)XMSE8M 
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WRITE(6,*)'MSELS' 

WRITE(6,*)XMSELS 

STOP 

END 
p     ****************  snRRnijTTMP  MSEBFT  **************************** 

C 

SUBROUTINE  MSEBET(BX,BY,AN,XMSEB) 

DIMENSION  BX(5000),BY(5000),SUM(5000) 

PLACE=0 

DO  901  1=1, AN 

SUM(I)=((BX(I)-1)*(BY(I)-1))**2 

PLACE=PLACE  +  SUM(I) 
901   CONTINUE 

XMSEB= PLACE/AN 

RETURN 

END 
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